We take the expression which appear in the inequalities of the domain to define new the variables,

({{u=(x^13)(y^2),}, {v=(x^32)(y^5)}})   ⟹  D≡({{{, ({{9<=u<=13,  }, {1<=v<=8,  }})}})

To compute the transformation that we need for the chain rule, we have to solve in the equations above for x and y,

Solve[{u==(x^13)(y^2),v==(x^32)(y^5)},{x,y}]

Therefore, the transformation is

T[u_,v_]:={((u^5)/(v^2)),((v^13)/(u^32))}

and the Jacobian determinant

Det[D[T[u,v],{{u,v},1}]]

The traduction of the function for which we have to compute the integral to the new variables is

2x+2y/.{x->((u^5)/(v^2)),y->((v^13)/(u^32))}

Finally the integral is

(∫_9^13)(∫_1^8)(((2 (u^5))/(v^2))+((2 (v^13))/(u^32)))((v^10)/(u^28))\[DifferentialD]v\[DifferentialD]u//N

We can represente the boundary of the domain consisting of the plane and the paraboloid indicated.

ContourPlot3D[{6x+2z==1,z==7(x^2)+7(y^2)},{x,-1,1},{y,-1,1},{z,0,3}]

Therefore, we can see that the volume indicated is the one over the paraboloid and under the plane (the other volume over the plane and inside the paraboloid is unbounded and with bolume ∞). Therefore, the domain is

D≡({{{, ({{6x+2z<=1,,  }, {z>=7(x^2)+7(y^2),  }})}})and the boundary     ∂D≡({{{, ({{6x+2z=1,,  }, {z=7(x^2)+7(y^2).,  }})}})

We can represent the domain by means of

RegionPlot3D[6x+2z<=1&&z>=7(x^2)+7(y^2),{x,-1,0.2},{y,-0.5,0.5},{z,0,2.5},PlotPoints->100,Boxed->False,Axes->True,AxesOrigin->{0,0,0}]

We can use cylindrical change of variable and then

({{x=ρ cos(θ)}, {y=ρ sin(θ)}, {z=z}, {(ρ^2) = (x^2) + (y^2)}})

If we translate to cylindrical coordinates of the domain and the boundary, we have

D≡({{{, ({{6ρ cos(θ)+2z<=1,,  }, {z>=7(ρ^2).,  }})}})       ∂D≡({{{, ({{6ρ cos(θ)+2z=1,,  }, {z=7(ρ^2).,  }})}})

For the limits in cylindrical coordinates we have that, as the domain is all arround the vertical axis the angle θ will move from 0 to 2π. For every angle θ, the radious ρ will move from 0 to the max distance to the vertical axis that we obtain for the intersection of the plane and the paraboloid and we can solve for ρ in terms of θ by means of the equations of the boundary:

Solve[{6ρ Cos[θ]+2z==1,z==7(ρ^2)},{ρ,z}]

As ρ is positive we take the only positive solution

ρ->(1/14) (-3 Cos[θ]+(Sqrt[14+9 ((Cos[θ])^2)]))

and therefore the radious ρ will move from 0 to (1/14) (-3 Cos[θ]+(Sqrt[14+9 ((Cos[θ])^2)])).

Once we know the angle and the radius, we can use the inequalities of the domain in cylindrical coordinates to determine that

D≡({{{, ({{6ρ cos(θ)+2z<=1,,  }, {z>=7(ρ^2).,  }})}})  ⟹  ({{{, ({{z<=(1/2)-3ρ cos(θ),,  }, {z>=7(ρ^2).,  }})}})

Therefore the limits for the integral in cyindrical coordinates are

(∫_0^(2π))(∫_0^((1/14) (-3 Cos[θ]+(Sqrt[14+9 ((Cos[θ])^2)]))))(∫_(7(ρ^2))^((1/2)-3ρ Cos[θ]))ρ\[DifferentialD]z\[DifferentialD]ρ\[DifferentialD]θ

N[%]

If we divide in a suitable way each inequality defining the domain we can recrite it as:

({{2(x^4)(y^6)<=z<=8(x^4)(y^6)}, {(x^3)(y^8)(z^9)<=1<=9(x^3)(y^8)(z^9)}, {4(x^8)(y^4)<=(z^8)<=13(x^8)(y^4)}})  ⟹  ({{2<=(z/((x^4)(y^6)))<=8}, {1<=(1/((x^3)(y^8)(z^9)))<=9}, {4<=((z^8)/((x^8)(y^4)))<=13}})

Therefore, if we take the new variables u, v, w defined by

({{u=(z/((x^4)(y^6))),}, {v=(1/((x^3)(y^8)(z^9))),}, {w=((z^8)/((x^8)(y^4))),}})

the limits for the domain would be

({{2<=u<=8,}, {1<=v<=9,}, {4<=w<=13.}})

Since we need the inverse transformation to apply the change rule for the integral, we have to solve in the equations defining u, v, w, the old variables x, y ,z. Thus,

Solve[{u==(z/((x^4)(y^6))),v==(1/((x^3)(y^8)(z^9))),w==((z^8)/((x^8)(y^4)))},{x,y,z}]

From all the possible solutions we have to select only the non complex and positive one which is:

{x->((u^(25/57))/((v^(11/57)) (w^(31/114)))),y->(((v^(2/19)) (w^(13/76)))/(u^(8/19))),z->((u^(13/57))/((v^(8/57)) (w^(7/114))))}

So we can define the transformation that we need from the simple domain in the variables u, v, w to the new ones:

T[u_,v_,w_]:={((u^(25/57))/((v^(11/57)) (w^(31/114)))),(((v^(2/19)) (w^(13/76)))/(u^(8/19))),((u^(13/57))/((v^(8/57)) (w^(7/114))))}

To apply the formula we also need the Jacobian determinant

Det[D[T[u,v,w],{{u,v,w},1}]]

and to translate the formula for which we have to compute the integral to the new variables:

(2x+(y^3))/.{x->((u^(25/57))/((v^(11/57)) (w^(31/114)))),y->(((v^(2/19)) (w^(13/76)))/(u^(8/19))),z->((u^(13/57))/((v^(8/57)) (w^(7/114))))}

Then, we finally have to compute the integral for new limits of the procut of the translated function with the Jacobian determinant:

(∫_2^8)(∫_1^9)(∫_4^13)(((2 (u^(25/57)))/((v^(11/57)) (w^(31/114))))+(((v^(6/19)) (w^(39/76)))/(u^(24/19))))(1/(228 (u^(43/57)) (v^(70/57)) (w^(265/228))))\[DifferentialD]w\[DifferentialD]v\[DifferentialD]u

N[%]